esort <- function(x, sortvar, ...) {
attach(x)
x <- x[with(x,order(sortvar,...)),]
return(x)
detach(x)
}

data <- read.csv("divergence.statistics.gimli.2.csv")


pdf("divergence_plots_b=0.02_d=0.016_nu=1e-04_t=7300.pdf",paper="letter",pointsize=8)
sw <- c(1e-06, 1e-07,1e-08)
pp <- c(0,0.25,0.5,0.75,1)
par(mfrow=c(3,5),mar=c(5,5,3,0.5))
for(i in 1:3){
for(ppi in 1:5){
dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==7300,]
ddmaxa0 <- dd[dd["max_min"]==1,]
ddmina0 <- dd[dd["max_min"]==-1,]
dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==7300,]
ddmaxa1 <- dd[dd["max_min"]==1,]
ddmina1 <- dd[dd["max_min"]==-1,]
plot(0,0,xlim=c(0,7300),ylim=c(0,0.3),pty="s",t="n",xaxs="r",yaxs="r",xlab="Time (years)",ylab="Divergence",las=1,xaxt="n")
axis(1,at=c(0,3650, 7300),lab=c(0,10,20))
title(paste("p=",pp[ppi]," mu=",sw[i],sep=""))
points(t(ddmaxa0["day"]),t(ddmaxa0["divergence"]),t="p",pch=1,col="black")
points(t(ddmina0["day"]),t(ddmina0["divergence"]),t="p",pch=1,col="red")
points(t(ddmaxa1["day"]),t(ddmaxa1["divergence"]),t="p",pch=19,col="black")
points(t(ddmina1["day"]),t(ddmina1["divergence"]),t="p",pch=19,col="red")
for(repi in 1:3){
dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==7300&data["rep"]==repi,]
dd <- dd[with(dd,order(day)),]
lines(t(dd["day"]),t(dd["divergence"]),lty=2)
}
for(repi in 1:3){
dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==7300&data["rep"]==repi,]
dd <- dd[with(dd,order(day)),]
lines(t(dd["day"]),t(dd["divergence"]),lty=1)
}
}
}
dev.off()


pdf("divergence_plots_b=0.02_d=0.016_nu=1e-04_t=14600.pdf",paper="letter",pointsize=8)
sw <- c(1e-06, 1e-07,1e-08)
pp <- c(0,0.25,0.5,0.75,1)
par(mfrow=c(3,5),mar=c(5,5,3,0.5))
for(i in 1:3){
for(ppi in 1:5){
dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==14600,]
ddmaxa0 <- dd[dd["max_min"]==1,]
ddmina0 <- dd[dd["max_min"]==-1,]
dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==14600,]
ddmaxa1 <- dd[dd["max_min"]==1,]
ddmina1 <- dd[dd["max_min"]==-1,]
plot(0,0,xlim=c(0, 14600),ylim=c(0,0.3),pty="s",t="n",xaxs="r",yaxs="r",xlab="Time (years)",ylab="Divergence",las=1,xaxt="n")
axis(1,at=c(0, 3650, 7300, 10950, 14600),lab=c(0,10,20,30,40))
title(paste("p=",pp[ppi]," mu=",sw[i],sep=""))
points(t(ddmaxa0["day"]),t(ddmaxa0["divergence"]),t="p",pch=1,col="black")
points(t(ddmina0["day"]),t(ddmina0["divergence"]),t="p",pch=1,col="red")
points(t(ddmaxa1["day"]),t(ddmaxa1["divergence"]),t="p",pch=19,col="black")
points(t(ddmina1["day"]),t(ddmina1["divergence"]),t="p",pch=19,col="red")
for(repi in 1:3){
dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==14600&data["rep"]==repi,]
dd <- dd[with(dd,order(day)),]
lines(t(dd["day"]),t(dd["divergence"]),lty=2)
}
for(repi in 1:3){
dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==14600&data["rep"]==repi,]
dd <- dd[with(dd,order(day)),]
lines(t(dd["day"]),t(dd["divergence"]),lty=1)
}
}
}
dev.off()




pdf("divergence_plots_b=0.02_d=0.005_nu=1e-04_t=14600.pdf",paper="letter",pointsize=8)
sw <- c(1e-06, 1e-07,1e-08)
pp <- c(0,0.25,0.5,0.75,1)
par(mfrow=c(3,5),mar=c(5,5,3,0.5))
for(i in 1:3){
for(ppi in 1:5){
dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.005&data["t"]==14600,]
ddmaxa0 <- dd[dd["max_min"]==1,]
ddmina0 <- dd[dd["max_min"]==-1,]
dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.005&data["t"]==14600,]
ddmaxa1 <- dd[dd["max_min"]==1,]
ddmina1 <- dd[dd["max_min"]==-1,]
plot(0,0,xlim=c(0, 14600),ylim=c(0,0.3),pty="s",t="n",xaxs="r",yaxs="r",xlab="Time (years)",ylab="Divergence",las=1,xaxt="n")
axis(1,at=c(0, 3650, 7300, 10950, 14600),lab=c(0,10,20,30,40))
title(paste("p=",pp[ppi]," mu=",sw[i],sep=""))
points(t(ddmaxa0["day"]),t(ddmaxa0["divergence"]),t="p",pch=1,col="black")
points(t(ddmina0["day"]),t(ddmina0["divergence"]),t="p",pch=1,col="red")
points(t(ddmaxa1["day"]),t(ddmaxa1["divergence"]),t="p",pch=19,col="black")
points(t(ddmina1["day"]),t(ddmina1["divergence"]),t="p",pch=19,col="red")
for(repi in 1:3){
dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.005&data["t"]==14600&data["rep"]==repi,]
dd <- dd[with(dd,order(day)),]
lines(t(dd["day"]),t(dd["divergence"]),lty=2)
}
for(repi in 1:3){
dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.005&data["t"]==14600&data["rep"]==repi,]
dd <- dd[with(dd,order(day)),]
lines(t(dd["day"]),t(dd["divergence"]),lty=1)
}
}
}
dev.off()










